Print["Mathematica is able to solve the Equations of Motion for a damped oscillator. Since this involves differential equations, we use DSolve. "]
Print["Another function I use is Timing. This is just to emphasize that Mathematica is faster than a human."]
Print["Below, Mathematica solves the most general case. Note that F = ma = m(Overscript[x, ..]) = - kx -c(OverDot[x]). We define (c/m)=2β, and (k/m)=((Subscript[ω, 0])^2)."]
DSolve[{x''[t]+2β x'[t]+((ω0)^2)x[t]==0},x[t],t][[1]]//Timing

Print["Below, Mathematica finds the coefficients C[1] and C[2] that satisfy initial conditions, e.g. position ((Subscript[x, 0])) and velocity ((Subscript[v, 0])). "]
DSolve[{x''[t]+2β x'[t]+((ω0)^2)x[t]==0,x[0]==x0,x'[0]==v0},x[t],t][[1]]//Timing
Print["This equation can be Simplified and FullSimplified, taking more and more computational power (see Timing). "]
DSolve[{x''[t]+2β x'[t]+((ω0)^2)x[t]==0,x[0]==x0,x'[0]==v0},x[t],t][[1]]//Simplify//Timing
DSolve[{x''[t]+2β x'[t]+((ω0)^2)x[t]==0,x[0]==x0,x'[0]==v0},x[t],t][[1]]//FullSimplify//Timing
Print["Note that the above equation gives (1/0) if (Subscript[ω, 0])=β, i.e. the critically damped oscillator."]
Print["The critically damped oscillator needs to be solved as a separate case, see below. "]
DSolve[{x''[t]+2ω0 x'[t]+((ω0)^2)x[t]==0,x[0]==x0,x'[0]==v0},x[t],t][[1]]//Timing

Manipulate[If[β1==ω0,Print["(Subscript[β, under]) = (Subscript[ω, 0]) caused an (1/0) error in the general oscillator equation. I reset the value of (Subscript[β, under]) to be (Subscript[ω, 0]) - 0.001. If β = (Subscript[ω, 0]), the expression for the critically damped oscillator should be used instead."];β1=ω0-0.001];If[β2==ω0,Print["(Subscript[β, over]) = (Subscript[ω, 0]) caused an (1/0) error in the general oscillator equation. I reset the value of (Subscript[β, over]) to be (Subscript[ω, 0]) + 0.001. If β = (Subscript[ω, 0]), the expression for the critically damped oscillator should be used instead."];β2=ω0+0.001];OscillatorManipul[t_,ω0_,β_,x0_,v0_]:=((E^(-t β)) (x0 (Sqrt[(β-ω0) (β+ω0)]) Cosh[t (Sqrt[(β^2)-((ω0)^2)])]+(v0+x0 β) Sinh[t (Sqrt[(β^2)-((ω0)^2)])]))/((Sqrt[(β-ω0) (β+ω0)]));
CriticallyDamped[t_,β_,x0_,v0_]:=(E^(-t β)) (x0+t (v0+x0 β));Plot[{OscillatorManipul[t,ω0,0,x0,v0],OscillatorManipul[t,ω0,β1,x0,v0],CriticallyDamped[t,ω0,x0,v0],OscillatorManipul[t,ω0,β2,x0,v0]},{t,0,tmax},PlotRange->All,PlotLegends->{"Undamped (β = 0)","Underdamped (β < (Subscript[ω, 0]))","Criticallydamped (β = (Subscript[ω, 0]))","Overdamped (β > (Subscript[ω, 0]))"},AxesLabel->{"t","x"}],{{tmax,12,"Plot time ((Subscript[t, max]))"},1,50,Appearance->"Labeled"},{{x0,1,"initial position ((Subscript[x, 0]))"},0,10,Appearance->"Labeled"},{{v0,0,"initial velocity ((Subscript[v, 0]))"},-10,10,Appearance->"Labeled"},{{ω0,1,"normal frequency ((Subscript[ω, 0]))"},1,10,Appearance->"Labeled"},{{β1,0.2,"underdamping friction ((Subscript[β, under]))"},0,ω0,Appearance->"Labeled"},{{β2,3,"overdamping friction ((Subscript[β, over]))"},ω0,20,Appearance->"Labeled"},SaveDefinitions->True,SynchronousUpdating->False,TrackedSymbols:>{tmax,x0,v0,ω0,β1,β2}]